%clear all
%close all

%profile on

units = 0;                          % selects the units, 0 for units of chi; 1 for MHz
Navg = 1000;                           % number of trajectories
dim = 8;
phi = 0.5*pi;                       %phase of LO

    
if units == 0
    
    eta = 0.5;
    
    Delta_q = 1600;     
    g = sqrt(Delta_q);              % corresponds to chi = 1
    chi = g^2/Delta_q;
    kappa1 = 2*chi;
    kappa2 = 2*chi;
    
    Delta_a = sqrt(3)*chi;
    Delta_b = -sqrt(3)*chi;
    
    gamma_d1 = 0/300;
    gamma_d2 = 0/300;
    gamma_d3 = 0/300;
    gamma_1 = 0.0;
    
  
    tmeas_on = 0/chi;
    tmeas_off = 20/(chi*eta);
    
    ton = tmeas_on + 6;
    toff = tmeas_off - 6;
    
    T = tmeas_off;                     % final time
    dt = 0.005;                     % time step
    N = T/dt;                       % number of time steps
 
    
    %epsilon0 = 2.0/sqrt(toff - ton); % This is to get the error probability = 0.025
    epsilon0 = 2.0/sqrt(16); % This is to get the error probability = 0.025
    
    sigma = 10;
    
    elseif units == 1
    
    Delta_q = 500e6;     
    g = sqrt(Delta_q*1e6);              % corresponds to chi = 1MHz
    chi = g^2/Delta_q;
    kappa1 = 2*chi;
    kappa2 = 2*chi;
    
    Delta_a = sqrt(3)*chi;
    Delta_b = -sqrt(3)*chi;
    
    gamma_d1 = 0.0;
    gamma_d2 = 0.0;
    gamma_d3 = 0.0;
    gamma_1 = 0.0;
    
    ton = 2/chi;
    toff = 18/chi;
  
    tmeas_on = 0/chi;
    tmeas_off = 18/chi;
    
    T = tmeas_off + 2;                     % final time
    dt = 0.005;                     % time step
    N = T/dt;                       % number of time steps
 
    
    epsilon0 = 1.5/sqrt(tmeas_off - tmeas_on); % This is to get the error probability = 0.025
    sigma = 10;
    
end

%% Allocate memory for the integrated signal and the populations

s_save = zeros(1,Navg);


overlap_p = zeros(1, Navg);
overlap_m = zeros(1, Navg);




%% Initilize the random number generator

%randn('state',sum(100*clock))  % changing seed
%randn('state',4563)  % changing seed This produced good res
%randn('state',523)  % changing seed
%randn('state',522343)  % changing seed



%% Generate the trajectories



%tic


% Initial state

initstate = 1/sqrt(dim)*ones(dim,1);
initMat = initstate*initstate';

init = reshape(transpose(initMat),dim^2,1);

ideal_p = 0.5*[1,0,0,1,0,1,1,0];
ideal_p = sparse(ideal_p'*ideal_p);
ideal_m = 0.5*[0,1,1,0,1,0,0,1];
ideal_m = sparse(ideal_m'*ideal_m);

rho_save = zeros(dim,dim,Navg);


parfor count = 1:Navg
%for count = 1:Navg
    
    disp(count)
    
    [rho,current,s,dWk,Am, int_field_000, int_field_011, alpha, beta] = get_trajectory(init, T, dt, gamma_1, gamma_d1, gamma_d2, gamma_d3, dim, kappa1, kappa2, epsilon0, chi, g, Delta_q, N, ton, toff, sigma, phi, tmeas_on, tmeas_off, Delta_a, Delta_b, eta );
    %[rho,current,s,dWk,Am, int_field_000, int_field_011] = get_trajectory_no_kron(init, T, dt, gamma_1, gamma_d1, gamma_d2, gamma_d3, dim, kappa1, kappa2, epsilon0, chi, g, Delta_q, N, ton, toff, sigma, phi, tmeas_on, tmeas_off, Delta_a, Delta_b, eta );

    s_save(count) = s;
    rho_qb = transpose(reshape(rho(:,end),dim,dim));
    rho_save(:,:,count) = rho_qb;
    
    if s < 0 
        overlap_m(count) = trace(rho_qb*ideal_m);
    end
    
    if s > 0
        overlap_p(count) = trace(rho_qb*ideal_p);
    end
    
    
    
%     for j = 1:N
%         
%         
%         rho_qb(:,:,j) = transpose(reshape(rho(:,j),dim,dim));
%         
%         Tr(j) = Tr(j) + trace(rho_qb(:,:,j));
%         purity(j) = purity(j) + trace(rho_qb(:,:,j)*rho_qb(:,:,j));
%         %overlap(j) = overlap(j) + trace(ideal_p*rho_qb(:,:,end));
%         
%         %rhotilde = YY*conj(rho_qb(:,:,j))*YY;
%         %eig_val = sort(eig(rho_qb(:,:,j)*rhotilde));
%         %C(j) = C(j) + max([0,real(sqrt(eig_val(4))) - real(sqrt(eig_val(3))) - real(sqrt(eig_val(2))) - real(sqrt(eig_val(1)))]);
%         
%         P_000(count,j) = rho_qb(1,1,j);
%         P_001(count,j) = rho_qb(2,2,j);
%         P_010(count,j) = rho_qb(3,3,j);
%         P_011(count,j) = rho_qb(4,4,j);
%         P_100(count,j) = rho_qb(5,5,j);
%         P_101(count,j) = rho_qb(6,6,j);
%         P_110(count,j) = rho_qb(7,7,j);
%         P_111(count,j) = rho_qb(8,8,j);
%         
%         C_000_001(count,j) = rho_qb(1,2,j);
%         C_000_011(count,j) = rho_qb(1,4,j);
%         C_000_111(count,j) = rho_qb(1,8,j);
%         
%     end
    
    %% Box-car filtering
    
%     count_th = 0;
%     
%     for sth = linspace(0,abs(s_p - s_m),Nth)
%         count_th = count_th + 1;
%         
%         if s(end) < s0 - sth %((s_p - sth < s(end)) && (s(end) < s_p + sth))
%             
%             rho_qb_p(:,:,count_th) = rho_qb_p(:,:,count_th) + rho_qb(:,:,end);
%             acc_p(count_th) = acc_p(count_th) + 1;
%         end
%         
%         if s(end) > s0 + sth %((s_m - sth < s(end)) && (s(end) < s_m + sth))
%             
%             rho_qb_m(:,:,count_th) = rho_qb_m(:,:,count_th) + rho_qb(:,:,end);
%             acc_m(count_th) = acc_m(count_th) + 1;
%         end
%         
%     end
    
end


%% looking at the density matrix for each value of threshold

% for i = 1:length(acc_p)
%
%     if acc_p(count_theta,i) > 0
%
%         rho_qb_p(:,:,i) = rho_qb_p(:,:,i)/acc_p(count_theta,i);
%         rho_qb_p_nofb(:,:,i) = rho_qb_p_nofb(:,:,i)/acc_p(count_theta,i);
%
%         purity_p(count_theta,i) = trace(rho_qb_p(:,:,i)*rho_qb_p(:,:,i));
%         purity_p_nofb(count_theta,i) = trace(rho_qb_p_nofb(:,:,i)*rho_qb_p_nofb(:,:,i));
%
%         overlap_p(count_theta,i) = trace(ideal_p*rho_qb_p(:,:,i));
%         overlap_p_nofb(count_theta,i) =  trace(ideal_p*rho_qb_p_nofb(:,:,i));
%
%         rhotilde = YY*conj(rho_qb_p(:,:,i))*YY;
%         eig_val = sort(eig(rho_qb_p(:,:,i)*rhotilde));
%         c_p(count_theta,i) = max([0,real(sqrt(eig_val(4))) - real(sqrt(eig_val(3))) - real(sqrt(eig_val(2))) - real(sqrt(eig_val(1)))]);
%
%         rhotilde = YY*conj(rho_qb_p_nofb(:,:,i))*YY;
%         eig_val = sort(eig(rho_qb_p_nofb(:,:,i)*rhotilde));
%         c_p_nofb(count_theta,i) = max([0,real(sqrt(eig_val(4))) - real(sqrt(eig_val(3))) - real(sqrt(eig_val(2))) - real(sqrt(eig_val(1)))]);
%
%
%
%     end
%
%     if acc_m(count_theta,i) > 0
%
%         rho_qb_m(:,:,i) = rho_qb_m(:,:,i)/acc_m(count_theta,i);
%         rho_qb_m_nofb(:,:,i) = rho_qb_m_nofb(:,:,i)/acc_m(count_theta,i);
%
%         purity_m(count_theta,i) = trace(rho_qb_m(:,:,i)*rho_qb_m(:,:,i));
%         purity_m_nofb(count_theta,i) = trace(rho_qb_m_nofb(:,:,i)*rho_qb_m_nofb(:,:,i));
%
%         overlap_m(count_theta,i) = trace(ideal_m*rho_qb_m(:,:,i));
%         overlap_m_nofb(count_theta,i) = trace(ideal_m*rho_qb_m_nofb(:,:,i));
%
%         rhotilde = YY*conj(rho_qb_m(:,:,i))*YY;
%         eig_val = sort(eig(rho_qb_m(:,:,i)*rhotilde));
%         c_m(count_theta,i) = max([0,real(sqrt(eig_val(4))) - real(sqrt(eig_val(3))) - real(sqrt(eig_val(2))) - real(sqrt(eig_val(1)))]);
%
%
%         rhotilde = YY*conj(rho_qb_m_nofb(:,:,i))*YY;
%         eig_val = sort(eig(rho_qb_m_nofb(:,:,i)*rhotilde));
%         c_m_nofb(count_theta,i) = max([0,real(sqrt(eig_val(4))) - real(sqrt(eig_val(3))) - real(sqrt(eig_val(2))) - real(sqrt(eig_val(1)))]);
%
%
%     end
%
% end


%% Total mean and variance

% P_000_mean = mean(P_000,1);
% P_001_mean = mean(P_001,1);
% P_010_mean = mean(P_010,1);
% P_011_mean = mean(P_011,1);
% P_100_mean = mean(P_100,1);
% P_101_mean = mean(P_101,1);
% P_110_mean = mean(P_110,1);
% P_111_mean = mean(P_111,1);
% 
% C_000_001_mean = mean(C_000_001,1);
% C_000_011_mean = mean(C_000_011,1);
% C_000_111_mean = mean(C_000_111,1);

%purity_save(count_theta) = purity(end)./Navg;
%overlap_save(count_theta) = overlap(end)./Navg;
%Tr_save(count_theta) = Tr(end)./Navg;
%C_save(count_theta) = C(end)./Navg;

% Tr_save = Tr_save/Navg;
%toc
%end

%save state_vs_eta_2.mat

